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Abstract. This paper is concerned with the optimal shape design of the newtonian 
^ I viscous incompressible fluids driven by the stationary nonhomogeneous Navier-Stokes 

" ' equations. We use three approaches to derive the structures of shape gradients for some 

given cost functionals. The first one is to use the Piola transformation and derive the 
. state derivative and its associated adjoint state; the second one is to use the differentia- 

I bility of a minimax formulation involving a Lagrangian functional with a function space 

(-H ! parametrization technique; the last one is to employ the differentiability of a minimax 

"j^ i formulation with a function space embedding technique. Finally we apply a gradient 

type algorithm to our problem and numerical examples show that our theory is useful 
for practical purpose and the proposed algorithm is feasible. 
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This paper deals with the optimal shape design for the stationary Navier-Stokes fiow. 
This problem is of great practical importance in the design and control of many in- 
dustrial devices such as aircraft wings, cars, turbines, boats, and so on. The control 
variable is the shape of the fluid domain, the object is to minimize some cost functionals 
that may be given by the designer, and flnally we can obtain the optimal shapes by 
numerical computation. 

Optimal shape design has received considerable attention already. Early works 
concerning on existence of solutions and differentiability of the quantity (such as, state, 
cost functional, etc.) with respect to shape deformation occupied most of the 1980s 
(see [2, 4, 14, 19]), the stabilization of structures using boundary variation technique 
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has been fully addressed in [4, 14, 19]. For the optimal shape design for Stokes flow, 
many people are contributed to it, such as O.Pironneau [13], J.Simon [16], ZM Gao 
et.al.[8, 9], and so on. 

In this paper, in order to derive the structures of shape gradients with respect to 
the shape of the variable domain for some given cost functionals in shape optimization 
problems for Navier-Stokes flow, we suggest the following three approaches: 

(i) use the Piola transformation and derive the state derivative with respect to the 
shape of the fluid domain and its associated adjoint state; 

(ii) utilize the differentiability of a minimax formulation involving a Lagrangian func- 
tional with a function space parametrization technique; 

(iii) employ the differentiability of a minimax formulation involving a Lagrangian 
functional with a function space embedding technique; 

In [7], we use the first approach to solve a shape optimization problem governed 
by a Robin problem, and in [9], we derive the expression of shape gradients for Stokes 
optimization problem by the first approach. In this paper, we use this approach to 
study the optimal shape design for Navier-Stokes flow with small regularity data. 

As we all known, many shape optimization problems can be expressed as a minimax 
of some suitable Lagrangian functional. Theorems on the differentiability of a saddle 
point (i.e., a minimax) of such Lagrangian functional with respect to a parameter pro- 
vides very powerful tools to obtain shape gradients by function space parametrization 
or function space embedding without the usual study of the state derivative approach. 

The function space parametrization technique and function space embedding tech- 
nique are advocated by M.C.Delfour and J.-P.Zolesio to solving poisson equation with 
Dirichlet and Nuemann condition (see[4]). In our paper [6, 8], we apply them to solve 
a Robin problem and a shape optimization problem for Stokes flow, respectively. How- 
ever, in this paper we extend them to study the optimal shape design for Navier-Stokes 
flow in despite of its lack of rigorous mathematical justification in case where the La- 
gragnian formulation is not convex. We shall show how this theorem allows, at least 
formally to bypass the study of the differentiability of the state and obtain the expres- 
sion of shape gradients with respect to the shape of the variable domain for some given 
cost functionals. 

We will find that the three approaches lead to the same expressions of the shape 
gradients for our given cost functionals. Hence, even if the last two approaches lacks 
from a rigorous mathematical framework, they allow more flexible computations which 
can be very useful for practical purpose. On the numerical point of view, we give the 
implementation of our problem in two dimensional case at the end of this paper, and 
the numerical results show that the last two approaches provide big efficiency for the 
shape optimization problem. 
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This paper is organized as follows. In section 2, we briefly recall the velocity method 
which is used for the characterization of the deformation of the shape of the domain 
and give the definitions of Eulerian derivative and shape derivative. We also give the 
description of the shape optimization problem for Navier-Stokes flow. 

In section 3, we prove the existence of the weak Piola material derivative, and give 
the description of the shape derivative. After that, we express the shape gradients of 
the typical cost functionals Ji{0,), (i = 1,2) by introducing the corresponding adjoint 
state systems. 

Section 4 is devoted to the computation of the shape gradient of the Lagrangian 
functional due to a minimax principle concerning the differentiability of the minimax 
formulation by function space parametrization technique and function space embedding 
technique. 

Finally in the last section, we give a gradient type algorithm with some numerical 
examples to prove that our theory could be very useful for the practical purpose and 
the proposed algorithm is feasible. 

2 Preliminaries and statement of the problem 

2.1 Elements of the velocity method 

Domains Q don't belong to a vector space and this requires the development of shape 
calculus to make sense of a "derivative" or a "gradient". To realize it, there are about 
three types of techniques: J.Hadamard [10] 's normal variation method, the perturbation 
of the identity method by J.Simon [15] and the velocity method (see J.Cea[2] and J.- 
P.Zolesio[4, 18]). We will use the velocity method which contains the others. In that 
purpose, we choose an open set D in with the boundary dD piecewise C^, and a 
velocity space V eE'' ■.= {V e C([0, e]; M^)) : V-uqd = on dD}, where e is a 

small positive real number and X'''(Z),R^) denotes the space of all /c— times continuous 
differentiable functions with compact support contained in . The velocity field 

V{t){x) = V{t,x), xeD, t>0 

belongs to V''{D,R^) for each t. It can generate transformations 

Tt{V)X = x{t, X), t>0, X £D 

through the following dynamical system 

(^{t,X) = V{t,xit)) 
\x{0,X)=X 

with the initial value X given. We denote the "transformed domain" Tt{V){0,) by 
nt{V) at t > 0, and also set Tt := Tt{T). 

There exists an interval I = [0,6), < 6 < e, and a one-to-one map Tt from D onto 
D such that 
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(i) To = I; 

(ii) {t,x) ^ Tt{x) belongs to C^{r,C^{D; D)) with Tt{dD) = dD; 

(iii) {t,x) ^ Tf^ix) belongs to C{I;C^{D- D)). 

Such transformation are well studied in [4]. 

Furthermore, for sufficiently small t > 0, the Jacobian Jt is strictly positive: 

Jt{x) := det|DTt(2;)| = detDrt(x) > 0, (2.2) 

where DTt [x) denotes the Jacobian matrix of the transformation Tt evaluated at a point 
X (z D associated with the velocity field V. We will also use the following notation: 
DTj~^(x) is the inverse of the matrix DTt(x) , *DTf~^{x) is the transpose of the matrix 
DTj~^(x). These quantities also satisfy the following lemma. 

Lemma 2.1 For any V G , DT^ and Jt are invertible. Moreover, DTt, DT^~^ are 
inCmO,e];C''-^{D;R^''^)), and Jt, J^"^ are in C^{[{),e]-C''-^{D-R)) 

Now let J(0) be a real valued functional associated with any regular domain we 
say that this functional has a Eulerian derivative at Vi in the direction V if the limit 

t\o t V > ; 

exists. 

Furthermore, if the map 

V ^ dJ{n;V) : E'' ^ R 

is linear and continuous, we say that J is shape differentiable at 0. In the distribu- 
tional sense we have 

dJ(r2;F) = (VJ, F)-pfc(5_iR]V)/xi5fe(D,R'V). (2.3) 

When J has a Eulerian derivative, we say that VJ is the shape gradient of J at 0,. 

Before closing this subsection, we introduce the following functional spaces which 
will be used throughout this paper: 

H{div,n) := {u € L'^{n)^ : div u = in u • n = on dn}, 
H^{dw,n) := {u € H^{n)^ : divu = in O}, 
H^{dw,n) := {u £ H\n)'^ : divti = in Q, u\dn = 0}. 
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2.2 Statement of the shape optimization problem 

Let be the fluid domain in M.'^ {N = 2 or 3), and the boundary F := dQ. The fluid 
is described by its velocity y and pressure p satisfying the stationary Navier-Stokes 
equations: 

—aAy + Dy ■ y + = / mQ 

divy = in^^ (2.4) 

\y = g on r 

where a stands for the inverse of the Reynolds number whenever the variables are 
appropriately nondimensionalized, / denotes the given body force per unit mass, and 
g is the given velocity at the boundary F. 

For the existence and uniqueness of the solution of the nonhomogeneous Navier- 
Stokes system (2.4), we have the following results (see [17]). 



Theorem 2.1 We suppose that O is of class . For 
g e ii't( div, M^) := \ g e 



H2 



t>N\ 



N 



dwg = 



(2.5) 
(2.6) 



there exists at least one y G H^{diY ,Q) and a distribution p G L'^{Q) on such that 
(2.4) holds. Moreover, if a is sufficiently large and g in — norm is sufficiently small, 
there exists a unique solution {y,p) € H^{dw ,Q) x Lq{Q) of (2.4). In addition, if Q, 
is of class C^, we have {y,p) € {H^{d\Y ^VL) r\ H'^iyi)^) x H^{n). 



We are interested in solving the following minimization problem 

min JifO) = - [ \y — vA"^ dx, 

or 

min J2(r2) = — |curly|^dx. 
Q£0 2 Jq 

An example of the admissible set O is: 

0:={ncR^ : |f| = 1}, 
where F is the domain inside the closed boundary F and |F| is its volume or area in 2D. 



(2.7) 
(2.8) 



3 State derivative approach 

In this section, we use the Piola transformation to bypass the divergence free condition 
and then derive a weak material derivative by the weak implicit function theorem. Then 
we will derive the structure of the shape gradients of the cost functionals by introducing 
the adjoint state equations associated with the corresponding cost functional. 
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3.1 Piola material derivative 

In order to deal with the nonhomogeneous Dirichlet boundary condition on F, we take 
y = y — g, where y satisfies the fohowing homogeneous Navier-Stokes system 

' -aAy + Dy ■ y + By ■ g + Bg ■ y + Vp = F in 
divy = inO (3.1) 

^ = on r 

with F := f + aAg + Bg g £ [^^(M^)]^. 

We say that the function y £ //g ( div , Q) is called a weak solution of problem (3.1) 
if it satisfies 

{e{y),w) = 0, weH^{div,n) (3.2) 

with 

{e{y), w) := / (aD^ : Bw + By ■ y ■ w + By ■ g ■ w + Bg ■ y ■ w — F ■ w) dx. (3.3) 
Jn 

As we all known, the divergence free condition coming from the fact that the fluid has 
an homogeneous density and evolves as an incompressible flow is difficult to impose on 
the mathematical and numerical point of view. Therefore in order to work with the 
divergence free condition, we need to introduce the following lemma (see [1]). 

Lemma 3.1 The Piola transform 

if^{{Jt)-^BTfip)oTf^ 

is an isomorphism. 

Now by the transformation Tt, we consider the solution y^ defined on 0^ of the 
perturbed weak formulation: 

/ {aBy :Bwt+By-y-Wt+By-g-Wt+Bg-y-Wt-F-Wt)dx = 0, Vwj G i/g ( div , Oj), 

(3.4) 

and introduce = "^^^{y^),w^ = "^^^{wt) defined on $7. 

We replace y,Wt by ^'^(j/*), ^'((lu*) in the weak formulation (3.4): 

[ [aB{^t{y')) : B{^t{w')) + B{^t{y')) ■ "^tiy") ■ ^^M) + B{^t{y')) ■ g ■ ^>t{w') 

+Bg ■ ^^tii/^) ■ ^^tiw^) - F ■ ^^tiw*)] dx = 0, Viu* eH^{dw,n). 

By the transformation Tf and the following identities: 

BiTf')oTt = BTf'; 
Bi^oTf') = {B<p.BTi')oT,-'; 
{Bg)oTt = B{goTt)-BT~\ 
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we use a back transport in and obtain the following weak formulation 

(e(t,y*),'u;*) = 0, y-w^ e H^{dw,Q) (3.5) 

with the notations 

{e{t,v),w) :=a D{B{t)v) : [B{B{t)w) ■ A{t)]dx 
Jn 

+ [ B{B{t)v)-v-{B{t)w)dx + [ [B{B{t)v)-BTf^]-{goTt)-{DTtw)dx 
Jn Jn 

+ / B{goTt)-v iB{t)w) dx- [ (Fo Tt) ■ (BTtw) dx, (3.6) 
Jn Jn 

and 

A{t) := Jt'DTf^*DTt-^; B{t)T := Jf^BTt ■ r. 
Now we are interested in the derivability of the mapping 

t^y' = ^;\yt) : [0, e] ^ H^oidiv , Q) 

where e > is sufficiently small and j/* G Hq ( div , Q) is the solution of the state 
equation 

{e{t,v),w) = 0, Vw e H^{dW,n). (3.7) 

In order to prove the differentiability of with respect to i in a neighborhood of t = 0, 
there maybe two approaches: 

(i) analysis of the differential quotient: lini(j/* — y)/t; 

(ii) derivation of the local differentiability of the solution y associated to the implicit 
equation (3.2). 

We use the second approach. However, we can not use the classical implicit theorem, 
since it requires strong differentiability results in H^^ for our case. Then we introduce 
the following weak implicit function theorem (see[18]). 

Theorem 3.1 Let X, Y' he two Banach spaces, I an open bounded set in M, and 
consider the map 

{t,x) ^ e(t,x) : I X X ^Y' 
If the following hypothesis hold: 

(i) t I— > {e{t,x),y) is continuously differentiahle for any y G Y and{t,x) i— > {dte{t, x),y) 
is continuous; 

(ii) there exists u £ X such that u € C^'^{I;X) and e(t,u(t)) = 0, Vt G /; 
(Hi) X I— > e{t,x) is differentiahle and {t,x) i— > dxe{t,x) is continuous; 

(iv) there exists tQ £ I such that dxe{t, x)\^t^^x{to)) isomorphism from X to Y' , 
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the mapping 

t ^ u{t) : I ^ X 

is dijferentiable at t = to for the weak topology in X and its weak derivative u{t) is the 
solution of 

{d^e{to, u{to)) • u{to),y) + {dte{to,u{to)),y) =0, Vy G Y. 

Now we state the main theorem of this subsection concerning on the differentiabihty 
of J/* with respect to t. 



Theorem 3.2 The weak Piola material derivative 

:=hm^ = Um5^L!M^ 

i->0 t t^O t 

exists and is characterized by the following weak formulation: 



I.e., 



{dve{0,v)\^^y-ij^,w) + {dte{0,y),w) =0, Vtu G div , 0), (3i 



a I T)y^ : Dw dx + Dy^ ■ y ■ w + Dy ■ y^ ■ w + Dy^ ■ g ■ w + Dg ■ y^ ■ w dx 

Jq Jn 

= a / By -.DwdwVdx + a / By : [-B{dwVw) +B(DVw) -DwBVjdx 
Jn Jn 

+ a Bw:[-B{diYVy) + B{BVy) -ByBV]dx 
Jn 

+ / [B(BVy - div Vy) y w + By y ■ {BVw - div Vw)] dx 
Jn 

+ I [B{BVy - div Vy) - ByBV] ■g-wdx + [ By- {BgV + *BVg) ■ w dx 
Jn Jn 

+ / [B{BgV) - Bg div V + *BVBg] ■ y ■ w dx 
Jn 

- [ [D(/ + aAg + Bg ■ g)V + *BV{f + aAg + Bg-g)]-w dx. (3.9) 
Jn 

where y is the solution of the weak formulation (3.2). 

Proof. In order to apply Theorem 3.1, we need to verify the fom' hypothesis of 
Theorem 3.1 for the mapping 

(t, v) ^ e(t, v) : [0, e] x H^{div , Q) ^ H^{div , ft)'. 

To begin with, since the flow map Tt G C\[0,£]]C^{D,R^)) and Lemma 2.1, the 
mapping 

t ^ {e{t,v),w) : [0,e] ^ R 

is for any v,w e HQ{div On the other hand, since F G [L^(M^)]^, the 

mapping t F o Tt is only weakly differentiable in the space [ff~^(M^)]^, thus the 



mapping t ^ e{t, v) is weakly differentiable. We denote by dte{t, v) its weak derivative. 
Since we have the foUowing three identities, 

-^DTt = {BV{t)oTt)m; (3.10) 

-^^Jt = {dwV{t))oTtJt; (3.11) 

A^foTt) = {Bf-Vit))oTt, (3.12) 
the weak derivative dte{t, v) can be expressed as fohows, 

{dte{t, v),w) = a [ B{B'{t)v) : [D(5(t)w) • A{t)] dx 
Jn 

+ a [ B{B{t)v) : [B{B'{t)w) ■ A{t) + B{B{t)w) ■ A'{t)] dx 
Jn 

+ / [B{B'{t)v) ■ V ■ iB{t)w) + B{B{t)v) ■ v ■ iB'{t)w)] dx 
Jn 

+ [ [B{B'it)v)BTf^ - B{B{t)v) ■ {BTf\BV{t) o Tt)BTi')] ■ (g o Tt) ■ (BTtw) dx 
Jn 

+ / [B{B{t)v)BTt-'] • {[{Bg • V{t)) o Tt] ■ {BTtw) + {g o Tt) ■ [{BV{t) o Tt)BTtw]}dx 
Jn 

+ / [Bi{Bg ■ Vit)) oTt)-v {Bit)w) + B{g oTt)-v- {B'{t)w)] dx 
Jn 

- I [*BVit)F + BFV{t)] o TtBTt w dx, (3.13) 
Jn 



where the notation 



B'{t)T := I [B{t)T] = [BV{t)oTt - {divV{t) oTt)l]B{t)T; 
A'{t) := |A(t) = [diyVit)oTt-BTi^BV{t)oTt]A{t)-*[BTi^BV{t)oTtA{t)]. 

It is easy to check that the mapping (t, v) i— s- dte{t, v) is weakly continuous from 
[0, e] X H^{div , Q) to H^{div , Q)'. 

We set t = 0, use Tt\t=o = I and V{t)\t=o = V, then obtain 

{dte{0,v),w) 

= a / Bv.BwdwVdx + a / Bv: [-B{divVw) + B{BVw) -BwBV]dx 
Jn Jn 

+ a Bw : [-D( div Vv) + B{BVv) - BvBV] dx 
Jn 

+ / [B{BVv - div Vv) -v-w + Bv-v- {BVw - div Vw)] dx 
Jn 

+ / [B{BVv - div Vv) - BvBV] ■g-wdx+ / Bv ■ {BgV + *BVg) ■ w dx 
Jn Jn 

+ / [B{BgV) - Bg div V + *BVBg] -y-wdx- [ {BFV + *BVF) ■ w dx. (3.14) 
Jn Jn 



To verify (ii), we follow the same steps described in R.Dziri[5] to find an identity satisfied 
by y^^ — y^^ and prove that the solution y^ € i7Q(div , 0) of the weak formulation 

{e{t,v),w) =Q, G i^oH div , f]) 

is Lipschitz with respect to t. 

It is easy to find that v i— > e(t, v) is differentiable, and the derivative of e(t, v) with 
respect to v in the direction 5v is 

{dve{t,v)-5v,w) = [ aB{B{t)6v) : [B{B{t)w)A{t)]dx 

+ [ [I){B{t)Sv) ■ V ■ {B{t)w) + D{B{t)v) ■ 6v ■ {B{t)w)] dx 
Jn 

+ / {[D{B{t)Sv)DTf^] ■ {g o Tt) ■ (DTtw) + D{g o Tt) • 5v ■ {B{t)w)} dx. 
Obviously, dve{t,v) is continuous, and when we take t = 0, 

{dve{0, v) ■ 6v, w) = a T)6v : Dw dx + D{6v) ■ v ■ w + Dv ■ 6v ■ w 

Jn Jn 

+ D{6v) ■ g ■ w + Dg ■ 6v ■ w dx. (3.15) 

Furthermore, the mapping 5v dve{0, v)-6v is an isomorphism from iJQ( div , $7) to its 
dual. Indeed, this result follows from the uniqueness and existence of the Navier-Stokes 
system, i.e., Theorem 2.1. 

Finally, all the hypothesis are satisfied by (3.6), we can apply Theorem 3.1 to (3.6) 
and then use (3.14) and (3.15) to obtain (3.9). □ 

3.2 Shape derivative 

In this subsection, we will characterize the shape derivative y', i.e., the derivative of 
the state y with respect to the shape of the domain. 

Theorem 3.3 Assume that $7 is of class C"^ , y G ^fg(div,Q) solves the weak formu- 
lation (3.2) in and y^ € -ffQ(div,r2() solves the perturbed weak formulation (3.4) in 
Of, then the shape derivative 

-/ 1- yt-y 

V := lim 

exists and is characterized as the solution of 

' -aAy' + T)y' ■ y + T)y ■ y' + Dy' ■ g + Dg ■ y' + Vp' = in Q 
< divj/' = inn (3.16) 

^y' = -{By-n)Vn on F 

Proof. We recall that y^ € i?Q(div,Of) satisfies the following weak formulation 

/ {ctDyf- : Dw + Dy^ ■ y^ - w + Dy^ ■ g ■ w + Dg ■ yf ■ w) dx — F ■ w dx = (3.17) 

Jfit Jnt 
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for any w G ifp ( div , Qt)- 

To begin with, we introduce the following Hadamard formula (see [4, 19]) 
d f f OF f 

— F{t,x)dx= (t,x)dx+ F{t,x)V -ntdTt, (3.18) 

for a sufficiently smooth functional F : [0, r] x — > R. 

Now we set a function (p € ^{Q)^ and div(^ = in Q. Obviously when t is 
sufficiently small, (f belongs to the sobolev space i^g ( div , ). Hence we can use 
(3.18) to differentiate (3.17) with w = (p, 

/ {aBy' -.Bff + Dy' - y ■ (f + By -y' -(f + Dy' - g -(f + Bg-y' ■(p)dx 
Jn 

+ j {aBy -.Bif + By ■ y ■ If + By ■ g ■ if + Bg ■ y ■ (f - F ■ (p)Vn ds = 0. 

Since (f has a compact support, the boundary integral vanishes. Using integration by 
parts for the first term in the distributed integral, we obtain 

/ i-aAy' + By' -y + By-y' + By' ■g + Bg-y')-<pdx = 0. (3.19) 
Jn 

Then there exists some distribution p' € L^(0) such that 

-aAy' + By' y + By y' + By' g + Bg y' = -Vp' 

in the distributional sense in 

Now we recall that for each t, ^^^(y^) belongs to the Sobolev space i?o(div,il), 
then we can deduce that its material derivative vanishes on the boundary T. Thus we 
obtain the shape derivative of y at the boundary, 

y' = -By • F, on r 

Since y|r = 0, we have Dj/|r = By ■ n*n, and then 

y' = -{By-n)Vn on T. 

□ 

Remark 3.1 Notice that in Theorem 3.3, the pressure p' is the shape derivative of the 
pressure pt which was defined on Of . 

The shape derivative y' of the solution y of the original Navier-Stokes system (2.4) is 
given by y' = y', then we obtain the following corollary by substituting y' = y' and 
y = y-g into (3.16). 

Corollary 3.1 The shape derivative y' of the solution y of (2.4) exists and satisfies 
the following system 

' -aAy' + By' ■y + By y' + Vp' = in n; 
< divy' = inn; (3.20) 

^y' = {B{g-y)-n)Vn onF. 

Moreover, we have y' G H^{div ,n). 
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3.3 Adjoint state system and gradients of the cost functionals 



This subsection is devoted to the computation of the shape gradients for the cost 
functionals Ji(ri) and J2i^) by the adjoint method. 

For the cost functional Ji{^) = Jq^^Iv — dx, we have 

Theorem 3.4 Let be of class and the velocity V G E^, the shape gradient VJi 
of the cost functional Ji (O) can be expressed as 



VJi 



^{y- Vd? + a(D(z/ -9)-n)- (Di? • n) 



(3.21) 



where the adjoint state v G ffQ(div,il) satisfies the following linear adjoint system 



—aAv — Dv ■ y + *Dy ■ v + Vg = y — y^, in ^ 
AiYV = 0, in 
= 0, on r. 



(3.22) 



Proof. Since Ji{^) is differentiable with respect to y, and the state y is shape 
differentiable with respect to t, i.e., the shape derivative y' exists, we obtain Eulerian 
derivative of Ji{^) with respect to t. 



d Ji(17; V)= l^{y-ya)-y'dx + l^^\y- yfVn ds 



(3.23) 



by Hadamard formula (3.18). 

By Green formula, we have the following identity 

/ [{-aAy' + By' ■y + By y' + Vp') ■ w - divy'vr] dx 
Jn 

= / [{—aAw — Bw ■ y + *By ■ w + Vvr) ■ y' — p' div w] dx 
Jn 

+ J {y' ■ w){y ■ n) ds + J {aBw ■ n — nn) ■ y' ds + J {pn — aBy'n)-wds. (3.24) 

Now we define {v, q) G H^{div , x L'^{n) to be the solution of (3.22), use (3.20) and 
set (io,7r) = {v,q) in (3.24) to obtain 



/ {y - Vd) ■ y' = - / ("Dt; -n- qn) -y' ds. 
Jn Jr 



(3.25) 



Since y' = {B{g — y) ■ n)Vn on the boundary T and divy' = in 17, we obtain the 
Eulerian derivative of Ji{0,) from (3.23), 



dJi{n;V) 



^ly- yd? + " {^{y -a)-^)- {Bv ■ n) 



Vn ds. 



(3.26) 



Since the mapping V ^ dJi{n;V) is linear and continuous, we get the expression 
(3.21) for the shape gradient VJi by (2.3). □ 
For another typical cost functional J2{S^) = §7^1 curlyp dx, we have the following 
theorem. 
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Theorem 3.5 Let $7 be of class and the velocity V € E^, the cost functional J2{^) 
possesses the shape gradient VJ2 which can be expressed as 



VJo = a 



^|curlyp + (D(y — g) ■ n) ■ (Dv ■ n — curly A n) 



(3.27) 



where the adjoint state v G {{^{div ,Q) satisfies the following linear adjoint system 

—aAv — Dv ■ y + *Dy • v + Vq = —aAy, in Q 

dwv = 0, inn (3.28) 
V = 0, on r. 

Proof. The proof is similar to that of Theorem 3.4. Using Hadamard formula (3.18) 
for the cost functional J2, we obtain the Eulerian derivative 

dJ2{n;V) = a / curly • curly' dx+ / -| curly] V„ ds. (3.29) 
Jn Jr 2 

Then, we define {v,q) G H^{div,n) x L'^{n) to be the solution of (3.28), use (3.20) 

and set {w,7r) = {v,q) in (3.24) to obtain 

a / Ay • y' dx = / a{J)v ■ n) ■ y' ds. 
Jn Jr 

Applying the following vectorial Green formula 



(3.30) 



{if ■ Aijj + curl (f ■ curl tp + div (p div ■0) dx 



I 

Jn 

= J {ip ■ { curl tp An) + (f ■ n div ■?/') ds 
for the vector functions y G H^{diY , Q) and y' G div ,0,), we obtain 

„ f crl y ^ crl d. + „ / A, ^ ^' d. = „ / ( curt , A „) ■ (3.31) 

Jn Jn Jr 

Combining (3.29), (3.30) with (3.31), we obtain the Eulerian derivative 



dJ2{n;V) 



a 



-|curlyp + (D(y — g) ■ n) ■ (Dv ■ n — curly A n) 



Finally we arrive at the expression (3.27) for the shape gradient VJ2. 



Vnds. (3.32) 
□ 



4 Function space parametrization and function space em- 
bedding 

In this section, we restrict our study to the minimization problem (2.7), and problem 
(2.8) follows similarly. In section 3, we have used the local differentiability of the state 
with respect to the shape of the fluid domain and the associated adjoint system to 
derive the shape gradient of the given cost functional. However, we do not need to 
analyze the differentiability of the state in many cases. In this section we derive the 
structure of the shape gradient for the cost functional Ji{0,) = ^ JqIu — Z/dPdx by 
function space parametrization and function space embedding techniques in order to 
bypass the study of the state derivative. 
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4.1 A saddle point formulation 

In this subsection, we shall describe how to build an appropriate Lagrange functional 
that takes account into the divergence condition and the nonhomogeneous Dirichlet 
boundary condition. 

We set / G [iJ^(M^)]^ and g € H^/'^{div then introduce a Lagrange multi- 

plier /X and a functional 

L{n,y,p,v,q,ii) = / [{aAy-I)y-y-Vp + f)-v+diYyq]dx+ {y-g)-tids (4.1) 
Jn Jr 

for e Y{n) X Q{Q), {v,q) G P{Q) x Q{n), and /i € H~^/^{r)^ with 

Now we're interested in the following saddle point problem 

inf sup L{Q,y,p,v,q,fj,) 

{y,p)eY{n)xQ{n) {v,q,fj,)eP{Q)xQ{Q)xH-''/^{r)'^ 

The solution is characterized by the following systems: 

(i) The state {y,p) is the solution of the problem 

—aAy + Dy • y + Vp = f iiiQ 

divy = inO (4.2) 

y = g on r 

(ii) The adjoint state {v,q) is the solution of the problem 

—aAv — Dv ■ y + *Dy ■ v + S/q = inQ 
divu = inO (4.3) 

V = on F; 

(iii) The multiplier: fj, = aDv n — qn, on L. 
Hence we obtain the following new functional, 

L{U,y,p,v,q) = / [{aAy — Dy-y — 'Vp+f)-v+dwyq]dx+ / {y—g)-{aDv-n — qn)ds. 
Jn Jt 

To get rid of the boundary integral, the following identities are derived by Green for- 
mula, 

/ {y-g) - (Di> ■n)ds= / [{y - g) ■ Av + D(y - g) : Du] dx; 
Jv Jn 

/ (y-g) - {qn)ds= / [div {y - g)q + {y - g) -Vqjdx. 
Jr Jn 
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Thus we introduce the new Lagrangian associated with (2.4) and the cost functional 

= kln\y-yd\^dx: 

G{n,y,p,v,q) = ^ \y - y^pdx + / [{aAy - Dy ■ y - Vp + f) ■ v + div yq] dx 
^ Jn Jn 

+ a / [{y -g) ■ Av + I){y - g) -.Dvjdx- / [ div (y - g)q + {y - g) ■ Vq] dx. 
Now the minimization problem (2.7) can be expressed as the fohowing form 

min inf sup G(Q,y,p,v,q). 

neo {y,p)eY{Q)xP{n) {v,q)€P{n)xQ{n) 

We can use the minimax framework to avoid the study of the state derivative with 
respect to the shape of the domain. The Karusch-Kuhn- Tucker (KKT) conditions wih 
furnish the shape gradient of the cost functional Ji{0,) by using the adjoint system. 
To begin with, we derive the formulation of the adjoint system which was satisfied by 

For {p,v,q) £ Q{Q) x P{^) x Q{0,), G{^l,y,p,v,q) is differentiable with respect to 
y £ Y{Q) and we get 

dyG{n,y,p,v,q) ■ 6y = / {y-ya)-6ydx+ / [aA{6y) - D{6y) ■ y -Dy ■ 6y] ■ v dx 

Jn Jn 

+ a [6y ■ Av + B{6y) -.Bvjdx - / 6y-Vqdx, y6y eY{n). 
Jn Jn 

Integrating by parts, we obtain 

dyG{^},y,p,v,q) ■ 5y = {aAv + Bv ■ y - *I)y ■ v - Vq + y - y^) ■ 6y dx. (4.4) 

Jn 

Similarly for {y,v,q) € y{^) x -P(O) x Q{Q), G{i},y,p,v,q) is differentiable with 
respect to p G Q{^), and we have 

dpG{Q.,y,p,v,q)-5p= / —V{5p)-vdx= I 5p diwdx, y5p&Q{0,). (4.5) 

Jn Jn 

Hence, (4.4) and (4.5) lead to the following linear adjoint system 

—aAv — Dv ■ y + *Dy ■ v + Vq = y — yd inQ. 
div?; = in 17 (4.6) 

V = Q on T; 

Given a velocity field V G and transformed domain Q.t '■= Tt{^), our main task 
of this section is to get the limit 

limMzM (4,7) 



with 



j(t) = inf sup G{nt,yt,pt,vt,qt), (4.8) 

{yt,pt)eY{nt)xQ{nt) {Vt,qt)eP(nt)xQ{nt) 
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where {y^,pt) satisfies 

-aAy^ + By^ ■y^ + Vpt = f in VLt 

div?/t = inOt (4.9) 

yt = 9 on Tt 

and {vt,qt) satisfies 

' -aAvt - Bvt ■ y^ + *Dy^ ■ Vt + Vqt = yt-yd m 
< dwvt = in fit (4.10) 

t>t = onTj; 

Unfortunately, the Sobolev space Y{flt), Qi^t), and P{flt) depend on the parameter 
t, so we need a theorem to differentiate a saddle point with respect to the parameter t, 
and there are two techniques to get rid of it: 

• Function space parametrization technique; 

• Function space embedding technique. 

4.2 Function space parametrization 

This subsection is devoted to the function space parametrization, which consists in 
transporting the different quantities (such as, a cost functional) defined on the variable 
domain back into the reference domain Q which does not depend on the perturbation 
parameter t. Thus we can use differential calculus since the functionals involved are 
defined in a fixed domain Q with respect to the parameter t. 

We parameterize the functions in H"^{Qt)'^ by elements of H"^{fl)'^ through the 
transformation : 



o T^-i : H"\n)'^ H'^ifttY, integer m > 0. 



where "o" denotes the composition of the two maps and d is the dimension of the 
function ^p. 

Note that since Tt and T^^ are diffeomorphisms, it transforms the reference domain 
O (respectively, the boundary F) into the new domain fit (respectively, the boundary 
Ft of fit). This parametrization can not change the value of the saddle point. We can 
rewrite (4.8) as 

j{t)= inf ^^^^ sup G{flt.yoTi\poT^\voT^\qoT^'). {A.ll) 

iy,p)'^Y{n)xQ{n) {v,q)£P{n)xQin) 

It amounts to introducing the new Lagrangian for {y,p,v,q) S Y{fl) x Q{fl) x P{fl) x 

G{t, y,p,v, q) = G{nt,y o Ti\po Ti\v o Ti\q o r,-^). 
The expression for G{t,y,p,v,q) is given by 

G{t,y,p,v,q) :=I,{t)+l2{t)+Ut)+h{t), (4.12) 
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where 



'dx: 



h{t):=l [ \yoTi'-y^ 

hit) := / [(aA(y o Tf') - D{y o r^^) • {y o Tf^) - V(p o T^^) + /] • o T, 

+ div {yoTi^){qoTi^)dx- 
hit) :=a [ [{y o Tf^ - g ■ A{v o Tf') + D(y o Tf' - g) : D(^ o r,"^)] dx; 

h{t) :=- [ [div (y o Tf' -g){qo Tf') + (y o Tf' - g) • V{q o T,"^)] dx, 
in* 

Now we introduce the theorem concerning on the differentiabihty of a saddle point (or 
a minimax). To begin with, some notations are given as fohows. 
Define a functional 

g : [0, r] X X X y ^ M 

with T > 0, and X, Y are the two topological spaces. 
For any t G [0, r], define 

g{t) = inf supC/(i,x,y) 

and the sets 

X{t) = {x* G X : g{t) = sn^y^y 0(.t, x*, y)} 
Y{t, x) = {y* G y : g{t, x, y') = supy^y G{t, x, y)] 

Similarly, we can define dual functionals 

h{t) = sup inf Q(t,x,y) 

and the corresponding sets 

Y{t) = {y* G y : h{t) = inf^ex Git, x, y*)} 
X{t, 2/) = {x* G X : g{t, x\y) = \nU^x Q{t, x, y)} 

Furthermore, we introduce the set of saddle points 

S{t) = {(x,y) G X X y : g{t) = g{t,x,y) = h{t)} 

Now we can introduce the following theorem (see [3] or page 427 of [4]): 

Theorem 4.1 Assume that the following hypothesis hold: 

(HI) S{t)^%, t G [0,t]; 

(H2) The partial derivative dtQ^t, x,y) exists in [0,r] for all 



(x,y) G 



teo,T] 



u 



X(0) X U Y{t) 
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(H3) There exists a topology Tx on X such that for any sequence : t.„ G [0, r]} with 
lim tn = 0, there exists x° G -'^(O) and a subsequence {inj.}, and for each k > 1, 

n yoo 

there exists Xn^ G such that 

(i) lim the Tx topology, 

n yoo 

\iminf dtg{t,xn„y) > dtg{0,x'',y), Vy G y(0); 

('iJ^J There exists a topology Ty onY such that for any sequence {tn '■ tn G [0, r]} with 
lim tn = 0, there exists y^ G Y{0) and a subsequence {tn,,,}, cind for each k > 1, 

n yoo 

there exists yn^. G Y{tn,,) such that 
(i) lim yn,, = y° in the Ty topology, 

n yoo 

(n) 

limsup9t^(t,x,y„J < 5*^(0, x, Vx G X(0). 

T/ien i/iere exists G X(0) x y(0) suc/i that 

dj(0) = lim»<a^ 

^ t\0 t 

= inf sup (9t^(0,x,?/) = 9(^(0, = sup inf dtg{0,x,y) (4.13) 

T/iis means that (x'^jy^) G ^(0) x y(0) is a saddle point of dtG{0,x,y). 

Following Theorem 4.1, we need to differentiate the perturbed Lagrange functional 
G{t, y,p, V, q). Since {y,p,v, q) G i^)^ x (Q) x (Q)^ x (n) provided that T is 
at less (see [17]), we can use Hadamard formula (3.18) to differentiate G{t, y,p, v, q) 
with respect to the parameter t > 0, 

dtGit,y,p,v,q) = I[{0)+I!^{0) + I'M+I'M, 

where 

/((O) := ^ {h{t)} = Jjy - y,) ■ (-Dy V) dx + \j^\y- yfVn ds; (4.14) 

1^(0) := 4 {h{t)] = [ {aAy -Byy-Vp + f)- {-BvV) dx 
ot t=o 

+ [ [aA{-ByV) - D(-DyF) y-By- (-ByV) + V(Vp ■V)]-v dx 
Jn 

+ I [dw{-ByV)q+ dwy{-Vq-V)]dx 
Jn 

+ j^[{aAy - By ■ y -Vp + f) ■ V + divyg]F„ds; (4.15) 
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4(0) := ^ {hit)} =a [ {{-ByV) ■ Av + [y - g) ■ A{-I)vV) 
+ B{-ByV) : Bv + B{y - g) : B{-BvV)} Ax 

+ a j[{y - g) ■ Av) + B{y - g) : Bv]Vn ds; (4.16) 

I'M ■■= jrAh{t)] =- / [d\w{-ByV)q+ diY{y-g){-Vq-V) 

+ i-DyV) ■Vq-{y-g)-V{Vq- V)] dx-J[ div (y - g)q +{y-g)- Vq]Vn ds. 

(4.17) 

Since (y,p) satisfies (2.4), divu = and v\y = 0, also by Green formula we can simplify 
(4.15) to 

1^(0) = f [aA{-ByV) - B[-ByV) y-By- {-ByV)] ■ v dx 
Jn 

+ / div{~ByV)qdx. (4.18) 
Jn 

By Green formula, y\r = g and v\t = 0, we can simplify (4.16) to 

/^(0)=a / {-ByV) ■ Av dx - a / A{-ByV) ■ v dx 
Jn Jn 

+ a j^[B{y-g):Bv]Vnds; (4.19) 

By Green formula, y\r = g and divy = divg = 0, (4.17) can be simplified to 

/^(O) = - [ [{-ByV) ■ Vq + div {-ByV)q] dx (4.20) 
Jn 

Adding (4.14), (4.18), (4.19) and (4.20) together, 

i=l 

= - [ [aAv • {ByV) - B{ByV) y-By {ByV) - {ByV) ■Vq + {y- y^) ■ {ByV)] dx 

+ ^ |^|y - yf + aB{y - g) : Bv"^ V-nds. (4.21) 

Since {v, q) are characterized by (4.6), we multiply the first equation of (4.6) by {ByV) 
and integrate over Q, then the distributional integral in (4.21) vanishes, finally we 
obtain the boundary expression for the Eulerian derivative 

dJi{n;V) = j^!^^\y-yf + aB{y-g):Bv^ V-nds. (4.22) 

Since y\r = 9 and i>|r = 0, we have B{y — g)\r = B{y — g) ■ n*n and Bv\y = Bv • n*n, 
thus we obtain an expression for the shape gradient 

VJi = l^ly - yA^ + «[D(l/ -9)-n]- {Bv • n)| n, (4.23) 

which is the same as (3.21) in the previous section. 
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4.3 Function space embedding 

In the previous subsection, we have used the technique of function space parametriza- 
tion in order to get the derivative oi j{t), i.e., 

= , , inf , , sup G{nt,yt,Pt,Vt,qt)- (4.24) 

with respect to the parameter t > 0. This subsection is devoted to a different method 
based on function space embedding technique. It means that the state and adjoint 
state are defined on a large enough domain D (called a hold- all [4]) which contains all 
the transformations {ilt '■ < t < e} of the reference domain for some small e > 0. 
For convenience, let D = M^. Use the function space embedding technique, 

j{t)= inf sup G{Qt,y,V,V,Q). (4.25) 



where the new Lagrangian 

G{Qt,y,v,v,Q) = l [ \y-yfdx+f [{a^y-i}y-y-vv+f)-v+d\vyQ]dx 

^ JUt Jnt 

+ a f [{y-g)-AV + Diy-g):DV]dx-[ [div {y - g)Q + {y - g) -VQjdx. 
Jut Jnt 

(4.26) 

Since G H\R^)^,ge H^^^{R^)^, and is of class C^, the solution {yt,Pt,Vt,qt) 

of (...) belongs to ii'3(j7^)iv ^ (n^) n Ll{nt)) x [H^nt]^ n H^{nt)^) x {H\nt)n 

Lg(J^i)) instead of Y{nt) x Q(Ot) x P(J7t) x Q{nt). Therefore, the sets 

X = Y = H^{R^)^ X H^-{R^), 

and the saddle points S{t) = X{t) x Y{t) are given by 

Xit) = {{y,V)eX: y\nt = yt, V\nt = Pt} (4.27) 
Y{t) = {(V, Q)eY: VU = ^t, Qb, = qt} (4.28) 

Using Theorem 4.1, we may make the conjecture that we can bypass the inf-sup and 
state 

dJ(0;F)= inf sup dtG{nt,y,V,V,Q)\t=o- (4.29) 

{y,v)ex{o) (V,Q)ey{0) 

Now we compute the partial derivative of the expression (4.26) by Hadamard formula 
(3.18), 

dtG{nt,y,v,v,Q)= [ [Wi{y,v) + W2{y,v,v,Q)]v -ntds {am) 

where 

W,{y,V):=^\y-yf + aB{y-g):I)V; 
m{y,V,V, Q) := {aAy -By-y-VV + f)-V + {y-g)- (aAV - VQ) + Qdivg. 
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and rit denotes the outward unit normal to the boundary Ft- 

We note that the expression (4.30) is a boundary integral on Tt which will not 
depend on {y,V) and (V, Q) outside of 0^, so the inf and the sup in (4.29) can be 
dropped, we then get 

dJiin- V) = jm{y, v) + W2{y,p, V, q)]Vn ds 

However, {y — g)\r = v\r = and divgi = imply W2{y,p,v, g) = on the boundary 
r. Finally we have 

dJiin; y) = ^ - yf + aB{y - g) : Di;| F„ ds 

As in the previous subsection, we also have the shape gradient of the functional J(r2), 

VJi = l^jy - yf + a[D(y -g)-n]- {Dv • n)| n, (4.31) 

which is the same as the expressions (3.21) and (4.23). 



5 Gradient algorithm and numerical simulation 

In this section, we will give a gradient type algorithm and some numerical examples in 
two dimensions to prove that our previous methods could be very useful and efficient 
for the numerical implementation of the shape optimization problems for Navier-Stokes 
flow. 



5.1 A gradient type algorithm 

In this subsection, we will describe a gradient type algorithm for the minimization of a 
cost function J{Vt). As we have just seen, the general form of its Eulerian derivative is 

dJ{^l■,V) = j VJ -Vds, 

where V J denotes the shape gradient of the cost functional J. Ignoring regularization, 
a descent direction is found by defining 

V = -hkVJ (5.1) 

and then we can update the shape 0, as 

nk = ii + hkV)n (5.2) 

where is a descent step at k-th iteration. 

There are also other choices for the definition of the descent direction. Since the 
gradient of the functional has necessarily less regularity than the parameter, an iterative 
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scheme like the method of descent deteriortates the regularity of the optimized param- 
eter. We need to project or smooth the variation into Hence, the method 
used in this paper is to change the scalar product with respect to which we compute a 
descent direction, for instance, H^{^1)'^. In this case, the descent direction is the unique 
element d G such that for every V G H^{Q)'^, 

[ Bd:BVdx= dJ{n;V). (5.3) 
Jn 

The computation of d can also be interpreted as a regularization of the shape gra- 
dient, and the choice of H^{fl)'^ as space of variations is more dictated by technical 
considerations rather than theoretical ones. 

The resulting algorithm can be summarized as follows: 

(1) Choose an initial shape r^oi 

(2) Compute the state system and adjoint state system, then we can evaluate the 
descent direction dk by using (5.3) with $7 = fi^; 

(3) Set ^k+i = (Id — hkdk) O^, where hj. is a small positive real number. 

The choice of the descent step hk is not an easy task. Too big, the algorithm is 
unstable; too small, the rate of convergence is insignificant. In order to refresh hk, we 
compare hk with hk-i- If {dk,dk~i)H^ is negative, we should reduce the step; on the 
other hand, if d^ and d^-i are very close, we increase the step. In addition, if reversed 
triangles are appeared when moving the mesh, we also need to reduce the step. 

In our algorithm, we do not choose any stopping criterion. A classical stopping 
criterion is to find that whether the shape gradients V J in some suitable norm is small 
enough. However, since we use the continuous shape gradients, it's hopeless for us to 
expect very small gradient norm because of numerical discretization errors. Instead, 
we fix the number of iterations. If it is too small, we can restart it with the previous 
final shape as the initial shape. 

5.2 Numerical examples 

To illustrate the theory, we want to solve the following minimization problem 

min^ 1 \y - yrfpdx (5.4) 



n 2 



n 



subject to 



—aAy + Dy ■ y + Vp = f inU 

dwy = inn (5.5) 

y = on r := Tout U Tin; 

The domain Q is an annuli, and its boundary has two parts: the outer boundary Tout 
is an unit circle which is fixed; the inner boundary Tin which is to be optimized. We 
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choose the body force / = (/i, /2) as follows: 
ayilbx"^ + 15y2 - 1) 



/i 



5(x2 + y2)3/2 



.^ + 4 1426 + 



31 



a;2 + y2 



753 



\/a;2 + y2 



1860 Vx2 + y2 



02/(15x2 + 15y2 - 1) 
5(x2 + y2)3/2 

- y 



2 1 / , 31 753 

x^ + —\ 1426 + ^ ^ + ^ 

775 I a;2 + 2/2 ^3.2 + 2 



1860^7x2 + y2 



The target velocity y^^ is determined by the data / and the target shape of the domain 

n. 

In this model problem, we have the following Eulerian derivative: 

dJ{n;V) = j^ !^^\y-yf + aBy:Bv^Vnds. 

We will solve this shape problem with two different target shapes: 
Case 1: A circle: Tin = {{x,y)\x'^ + 2/^ = 0.22}. 
Case 2: An ellipse: T,^ = {{x,y)\^ + ^ = ^}. 

Our numerical solutions are obtained under FreeFem++ [12] and we run the program 
on a home PC. 

In Case 1, we choose the initial shape to be elliptic: |(x,2/) ^ + ^ = j^|i ^-^id the 
initial mesh is shown in Figure 5.1. 

In Case 2, we take the initial shape to be a circle whose center is at origin with 
radius 0.6, and the initial mesh is shown in Figure 5.2. 





Figure 5.1: Initial mesh in 
Case 1 with 226 nodes. 



Figure 5.2: Initial mesh in 
Case 2 with 161 nodes. 



In Case 1, Figure 5.3 — Figure 5.6 give the comparison between the target shape 
with iterated shape for the viscosity coefficients a = 0.1,0.01 and 0.001, respectively. 
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For a = 0.1 and a = 0.01, we choose the initial step /i = 20 and in Figure 5.3 
and Figure 5.4, we give the final shape at iteration 10 with CPU times. However for 
a = 0.001, one can not obtain a good result when h = 20 (see Figure 5.5). Thus we 
should reduce the initial descent step, and then in Figure 5.6, we give the final shape 
at iteration 40 with the initial step h = 5. By comparison with Figure 5.5, we find that 
though we need much more CPU time for h = 5, but it have a nicer reconstruction. 

Figure 5.7 represents the fast convergence of the cost functional for the various 
viscosities in Case 1. 




5 10 15 20 25 30 35 40 

Iteration 



Figure 5.7: Convergence history in Case 1 for a = 0.1,0.01 and 0.001. 

In Case 2, Figure 5.8 and Figure 5.9 show the comparison between the target shape 
with iterated shape at iteration 15 for the viscosity a = 0.1 and a = 0.01, respectively. 
At this time, we choose the initial step h = 20. We also give the CPU run times at 
the 15 iterations for a = 0.1 and a = 0.01. Unfortunately, we can not get a good 
reconstruction for a = 0.001 in this case. Figure 5.10 gives the convergence history of 
the cost functional J{Q) for a = 0.1,0.01. 

Finally, we can conclude that the proposed gradient type algorithm is an efficient one 
in both of our test cases. Unfortunately for large Reynold numbers, we can not obtain 
the nice results quickly. Hence further research is necessary on efficient implementations 
for very large Reynold numbers and real problems in the industry. 
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Figure 5.8: Case 2: a = 0.1, CPU time: 196.609 s. 
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Figure 5.9: Case 2: a = 0.01, CPU time: 207.469 s. 
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Figure 5.10: Convergence history in Case 2 for a = 0.1,0.01. 
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